library(readxl)
library(DT)
library(data.table)
library(dplyr)
library(ggplot2)
library(plotly)
library(cowplot)
library(ggrepel)
library(curl)
library(biomaRt)
library(sqldf)
# Ensembl LD API
library(httr)
library(jsonlite)
library(xml2)
library(gaston)
library(RCurl)
# *** susieR ****
# library(knitrBootstrap) #install_github('jimhester/knitrBootstrap')
library(susieR) # devtools::install_github("stephenslab/susieR")
# *** finemapr ****
## finemapr contains: finemap, CAVIAR, and PAINTOR
# library(finemapr) # devtools::install_github("variani/finemapr")
# library(roxygen2) #roxygenize()
# *** locuscomparer ****
# https://github.com/boxiangliu/locuscomparer
# library(locuscomparer); #devtools::install_github("boxiangliu/locuscomparer")
# thm <- knitr::knit_theme$get("bipolar")
# knitr::knit_theme$set(thm)
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] susieR_0.6.2.0411 RCurl_1.95-4.11 bitops_1.0-6
## [4] gaston_1.5.4 RcppParallel_4.4.2 Rcpp_1.0.0
## [7] xml2_1.2.0 jsonlite_1.6 httr_1.4.0
## [10] sqldf_0.4-11 RSQLite_2.1.1 gsubfn_0.7
## [13] proto_1.0.0 biomaRt_2.38.0 curl_3.3
## [16] ggrepel_0.8.0 cowplot_0.9.4 plotly_4.8.0
## [19] ggplot2_3.1.0 dplyr_0.8.0.1 data.table_1.12.0
## [22] DT_0.5.1 readxl_1.3.0
##
## loaded via a namespace (and not attached):
## [1] lattice_0.20-38 tidyr_0.8.2 prettyunits_1.0.2
## [4] assertthat_0.2.0 digest_0.6.18 R6_2.4.0
## [7] cellranger_1.1.0 plyr_1.8.4 chron_2.3-53
## [10] stats4_3.5.1 evaluate_0.13 pillar_1.3.1
## [13] rlang_0.3.1 progress_1.2.0 lazyeval_0.2.1
## [16] blob_1.1.1 S4Vectors_0.20.1 Matrix_1.2-15
## [19] rmarkdown_1.11 stringr_1.4.0 htmlwidgets_1.3
## [22] bit_1.1-14 munsell_0.5.0 compiler_3.5.1
## [25] xfun_0.5 pkgconfig_2.0.2 BiocGenerics_0.28.0
## [28] htmltools_0.3.6 tcltk_3.5.1 tidyselect_0.2.5
## [31] expm_0.999-3 tibble_2.0.1 IRanges_2.16.0
## [34] XML_3.98-1.17 viridisLite_0.3.0 crayon_1.3.4
## [37] withr_2.1.2 grid_3.5.1 gtable_0.2.0
## [40] DBI_1.0.0 magrittr_1.5 scales_1.0.0
## [43] stringi_1.3.1 tools_3.5.1 bit64_0.9-7
## [46] Biobase_2.42.0 glue_1.3.0 purrr_0.3.0
## [49] hms_0.4.2 parallel_3.5.1 yaml_2.2.0
## [52] AnnotationDbi_1.44.0 colorspace_1.4-0 memoise_1.1.0
## [55] knitr_1.21
print(paste("susieR ", packageVersion("susieR")))## [1] "susieR 0.6.2.411"
createDT <- function(DF, caption="", scrollY=400){
data <- DT::datatable(DF, caption=caption,
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = scrollY, scrollX=T, scrollCollapse = T, paging = F,
columnDefs = list(list(className = 'dt-center', targets = "_all"))
)
)
return(data)
}
createDT_html <- function(DF, caption="", scrollY=400){
print( htmltools::tagList( createDT(DF, caption, scrollY)) )
}
# https://stackoverflow.com/questions/42866055/rshiny-download-button-within-rmarkdown
# downloadHandler(filename = function() {
# return(paste('Example', input$SS, '.csv', sep=''))
#
# }, content = function(file) {
# write.csv(RandomSample(), file)
# })Use bioMart to get TSS positions for each gene
get_TSS_position <- function(gene){
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
# att <- listAttributes(mart)
# grep("transcription", att$name, value=TRUE)
TSS <- getBM(mart=mart,
attributes=c("hgnc_symbol","transcription_start_site", "version"),
filters="hgnc_symbol", values=gene)
} For each gene, get the position of top SNP (the one with the greatest effect size/Beta)
import_sig_GWAS <- function(file_path, caption="", sheet = 1,
chrom_col="CHR", position_col="POS", snp_col="SNP",
pval_col="P", effect_col="Effect", gene_col="Gene"){
## Only the significant subset of results
SumStats_sig <- read_excel(path = file_path, sheet = sheet)[,c(chrom_col, position_col, snp_col,
pval_col, effect_col, gene_col )]
# Standardize names
colnames(SumStats_sig) <- c("CHR","POS","SNP","P","Effect","Gene")
top_SNPs <- SumStats_sig %>% group_by(Gene) %>% top_n(-1, "P") #`Beta, all studies`
top_SNPs <- cbind(Coord=paste(top_SNPs$CHR, top_SNPs$POS, sep=":"),
top_SNPs)
createDT(SumStats_sig, caption)
return(list(top_SNPs, SumStats_sig))
}
# list[top_SNPs, SumStats_sig] <- import_sig_GWAS(
# file_path = "Data/Parkinsons/Nalls2018_S2_SummaryStats.xlsx",
# sheet="Data",
# chrom_col = "CHR", position_col = "BP", snp_col="SNP",
# pval_col="P, all studies", effect_col="Beta, all studies", gene_col="Nearest Gene",
# caption= "Nalls et al. (2018) PD GWAS Summary Stats")
# list[top_SNPs, SumStats_sig] <- import_sig_GWAS(
# file_path = "Data/Alzheimers/Posthuma/AD_target_SNP.xlsx",
# sheet = 3,
# chrom_col = "Chr", position_col = "bp", snp_col="SNP",
# pval_col="P", effect_col="Z", gene_col="Gene",
# caption= "Posthuma et al. (2018) AD GWAS Summary Stats")
# Nall (2018) QTL Stats
# Nalls_QTL <- read_excel("Parkinsons_Data/Nalls2018_S4_QTL-MR.xlsx")
# createDT(Nalls_QTL, caption = "Nalls (2018) Table S4. Expanded summary statistics for QTL Mendelian randomization") Get all genes surrounding the index SNP (default is 500kb upstream + 500kb downstream)
# 1000000 bp
get_flanking_SNPs <- function(gene, top_SNPs, bp_distance=500000, filePath,
chrom_col="CHR", position_col="POS", snp_col="SNP",
pval_col="P", effect_col="Effect", stderr_col="StdErr"){
read.delim(filePath, nrows = 2)
if(gene %in% top_SNPs$Gene==F){
cat("\n ----- Could not find gene '",gene,"' in topSNPs dataframe. Try a different gene name. ----- \n" )
}else{
topSNP_sub <- top_SNPs[top_SNPs$Gene==gene,]
minPos <- topSNP_sub$POS - bp_distance
maxPos <- topSNP_sub$POS + bp_distance
query <- sqldf::read.csv.sql(filePath, sep="\t", stringsAsFactors=F,
sql = paste('select * from file where',chrom_col,'=',topSNP_sub$CHR,
'AND',position_col,'BETWEEN', minPos, 'AND', maxPos))
geneSubset <- query %>% dplyr::rename(CHR=chrom_col,POS=position_col, SNP=snp_col,
P=pval_col, Effect=effect_col, StdErr=stderr_col) %>%
mutate(Effect=as.numeric(Effect), StdErr=as.numeric(StdErr), P=as.numeric(P))
## Remove SNPs with NAs in stats
geneSubset <- geneSubset[complete.cases(geneSubset),]
return(geneSubset)
# Close the connection
base::close(query)
unlink(filePath)
}
}
# flankingSNPs <- get_flanking_SNPs("LRRK2", top_SNPs,
# filePath = "Data/Parkinsons/META.PD.NALLS2014.PRS.tsv",
# snp_col = "MarkerName", pval_col = "P.value")
# flankingSNPs <- get_flanking_SNPs("CLU/PTK2B", top_SNPs,
# filePath = "Data/Alzheimers/Posthuma/phase3.beta.se.hrc.txt",
# effect_col = "BETA", stderr_col = "SE", position_col = "BP")gaston_LD <- function(flankingSNPs, reference="1KG_Phase1", superpopulation="EUR"){
# Download portion of vcf from 1KG website
region <- paste(flankingSNPs$CHR[1],":",min(flankingSNPs$POS),"-",max(flankingSNPs$POS), sep="")
chrom <- flankingSNPs$CHR[1]
if(reference=="1KG_Phase3"){
cat("LD Reference Panel = 1KG_Phase3 \n")
vcf_URL <- paste("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr",chrom,
".phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz",sep="")
if(!exists("popDat_1KGphase3")){
popDat <- popDat_1KGphase3 <- read.delim("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel", header = T, row.names = NULL)
colnames(popDat_1KGphase3) <- c("sample","population","superpopulation","gender")
} else{popDat <- popDat_1KGphase3}
} else if (reference=="1KG_Phase1") {
cat("LD Reference Panel = 1KG_Phase1 \n")
vcf_URL <- paste("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/ALL.chr",chrom,
".phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz", sep="")
# Download individual-level metadata
if(!exists("popDat_1KGphase1")){
popDat <- popDat_1KGphase1 <- read.delim("ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/phase1_integrated_calls.20101123.ALL.panel", header = F, row.names = NULL)
colnames(popDat) <- c("sample","population","superpopulation","seq_platform1","seq_platform2")
} else { popDat <- popDat_1KGphase1 }
}
# library(Rsamtools); #BiocManager::install("Rsamtools")
system(paste("tabix -fh ",vcf_URL ,region, "> subset.vcf"))
vcf_name <- paste(basename(vcf_URL), ".tbi",sep="")
# Import w/ gaston and further subset
bed.file <- read.vcf("subset.vcf")
## Subset rsIDs
bed <- select.snps(bed.file, id %in% flankingSNPs$SNP)
# Subset Individuals
selectedInds <- subset(popDat, superpopulation %in% superpopulation)
bed <- select.inds(bed, id %in% selectedInds$sample)
# Cleanup extra files
rm(bed.file)
file.remove(vcf_name)
file.remove("subset.vcf")
# Calculate pairwise LD for all SNP combinations
ld.x <- gaston::LD(bed, lim = c(1,ncol(bed)))
# LD plot
limit <- 20
LD.plot( ld.x[1:limit,1:limit], snp.positions = bed@snps$pos[1:limit])
# Double check subsetting
# LD_matrix <- ld.x[row.names(ld.x) %in% flankingSNPs$SNP, colnames(ld.x) %in% flankingSNPs$SNP]
LD_matrix <- ld.x
LD_matrix[!is.finite(LD_matrix)] <- 0
return(LD_matrix)
}
# LD_matrix <- gaston_LD(flankingSNPs)
# LD_matrix <- gaston_LD(flankingSNPs)susie_on_gene <- function(gene, top_SNPs,
bp_distance=500000, filePath, num_causal=1,
chrom_col="CHR", position_col="POS", snp_col="SNP",
pval_col="P", effect_col="Effect", stderr_col="StdErr",
LD_reference="1KG_Phase1", superpopulation="EUR"){
cat("\n + Extracting SNPs flanking lead SNP... \n")
flankingSNPs <- get_flanking_SNPs(gene, top_SNPs, bp_distance=bp_distance, filePath=filePath,
chrom_col=chrom_col, position_col=position_col, snp_col=snp_col,
pval_col=pval_col, effect_col=effect_col, stderr_col=stderr_col)
### Get LD matrix
cat("\n + Creating LD matrix... \n")
LD_matrix <- gaston_LD(flankingSNPs, LD_reference, superpopulation)
## Turn LD matrix into positive semi-definite matrix
# LD_matrix2 <- ifelse(matrixcalc::is.positive.semi.definite(LD_matrix),
# LDmatrix,
# Matrix::nearPD(LD_matrix)$mat %>% as.matrix() )
## Subset summary stats to only include SNPs found in query
geneSubset <- flankingSNPs %>% subset(SNP %in% unique(row.names(LD_matrix), colnames(LD_matrix) ) )
LD_matrix <- LD_matrix[ geneSubset$SNP, geneSubset$SNP]
b <- geneSubset$Effect
se <- geneSubset$StdErr
# Run Susie
cat("\n + Fine mapping with SusieR... \n")
fitted_bhat <- susie_bhat(bhat = b, shat = se,
R = LD_matrix,
n = nrow(LD_matrix),
L = num_causal, # 1
scaled_prior_variance = 0.1,
estimate_residual_variance = TRUE, # TRUE
estimate_prior_variance = FALSE, # FALSE
verbose = T,
# var_y = var(b),
standardize = TRUE
)
# Format results
geneSubset$Coord <- paste(geneSubset$CHR, geneSubset$POS, sep=":")
susieDF <- data.frame(SNP=names(fitted_bhat$X_column_scale_factors), PIP=fitted_bhat$pip) %>%
base::merge(subset(geneSubset, select=c("CHR","POS","SNP","Effect","P","Coord")), by="SNP") %>%
mutate(POS=as.numeric(POS))
return(susieDF)
}
# susieDF <- susie_on_gene("LRRK2", top_SNPs,
# filePath = "Data/Parkinsons/META.PD.NALLS2014.PRS.tsv",
# snp_col = "MarkerName", pval_col = "P.value")
# susieDF <- susie_on_gene(gene="CLU/PTK2B", top_SNPs,
# filePath="Data/Alzheimers/Posthuma/phase3.beta.se.hrc.txt",
# effect_col = "BETA", stderr_col = "SE", position_col = "BP")before_after_plots <- function(gene, susieDF, topVariants=3){
roundBreaks <- seq(plyr::round_any(min(susieDF$POS),10000), max(susieDF$POS),250000)
## Before fine-mapping
geneSubset <- susieDF %>% arrange(desc(abs(Effect)))
labelSNPs_Effect <- geneSubset[1:topVariants,]
yLimits <- c(min(-log(susieDF$Effect)), max(-log(susieDF$Effect))+1)
before_plot <- ggplot(geneSubset, aes(x=POS, y=Effect, label=SNP, color= -log(P) )) +
ylim(yLimits) +
geom_hline(yintercept=0,alpha=.5, linetype=1, size=.5) +
geom_point(alpha=.5) +
geom_point(data=labelSNPs_Effect[1,], pch=21, fill=NA, size=4, colour="red", stroke=1) +
geom_segment(aes(xend=POS, yend=0, color= -log(P)), alpha=.5) +
geom_text_repel(data=labelSNPs_Effect, aes(label=SNP), segment.alpha = .2, nudge_x = .5) +
labs(title=paste(gene," (",length(susieDF$PIP)," variants)","\nBefore Fine Mapping",sep=""), y="Beta", x="Position",
color="GWAS\n-log(p-value)") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = roundBreaks)
## After fine-mapping
susieDF <- susieDF %>% arrange(desc(PIP))
labelSNPs_PIP <- susieDF[1:topVariants,]
yLimits <- c(min(susieDF$PIP), max(susieDF$PIP)+.1)
after_plot <- ggplot(susieDF, aes(x=POS, y=PIP, label=SNP, color= -log(P) )) +
ylim(yLimits) +
geom_hline(yintercept=0,alpha=.5, linetype=1, size=.5) +
geom_point(alpha=.5) +
geom_point(data=susieDF[susieDF$PIP == max(susieDF$PIP),],
pch=21, fill=NA, size=4, colour="green", stroke=1) +
geom_segment(aes(xend=POS, yend=yLimits[1], color= -log(P))) +
geom_text_repel(data=labelSNPs_PIP, aes(label=SNP), segment.alpha = .2, nudge_x = .5) +
labs(title=paste(gene," (",length(susieDF$PIP)," variants)","\nAfter Fine Mapping",sep=""), y="PIP", x="Position",
color="GWAS\n-log(p-value)") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(breaks = roundBreaks)
plot_grid(before_plot, after_plot, nrow = 2) %>% print()
# susie_plot(fitted_bhat, y="PIP", b=b, add_bar = T)
createDT_html(susieDF, paste("susieR Results: ", gene), scrollY = 200)
}
# before_after_plots(gene = "LRRK2", susieDF, topVariants = 3)
# before_after_plots(gene = "CLU/PTK2B", susieDF, topVariants = 3)before_after_consensus <- function(gene, SumStats_sig, susieDF, max_SNPs=10){
# Get top SNPs from Nalls et all fine mapping
SS_dat <- subset(SumStats_sig, Gene==gene) %>%
arrange(P) %>% mutate(SNP=as.character(SNP))
SumStats_SNPs <-if(max_SNPs > length(SS_dat$SNP)){ SS_dat$SNP}else{SS_dat$SNP[1:max_SNPs]}
# Get top SNPs from susieR fine mapping
susieR_dat <- subset(susieDF, PIP!=0) %>%
arrange(desc(PIP)) %>% mutate(SNP=as.character(SNP))
susieR_SNPs <-if(max_SNPs > length(susieR_dat$SNP)){ susieR_dat$SNP}else{susieR_dat$SNP[1:max_SNPs]}
# Calculate percent
overlap <- length(intersect(SumStats_SNPs, susieR_SNPs))
percentOverlap <- overlap / length(susieR_SNPs) * 100
cat("\n",overlap," / ",length(susieR_SNPs), " (",round(percentOverlap,2),
"%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.","\n",sep="")
}
# before_after_consensus(SumStats_sig, susieDF, max_SNPs=10)#Nalls_SS %>% group_by(`Nearest Gene`) %>% tally() %>% subset(n>2)
finemap_geneList <- function(top_SNPs, geneList, filePath,
bp_distance=500000, num_causal=1,
chrom_col="CHR", position_col="POS", snp_col="SNP",
pval_col="P", effect_col="Effect", stderr_col="StdErr",
LD_reference="1KG_Phase1", superpopulation="EUR",
topVariants=3){
fineMapped_topSNPs <- data.table()
for (gene in geneList){
cat("\n")
cat(paste("###", gene, "\n"))
susieDF <- susie_on_gene(gene=gene, top_SNPs=top_SNPs, num_causal = 1,
filePath=filePath, bp_distance=bp_distance,
chrom_col=chrom_col, position_col=position_col, snp_col=snp_col,
pval_col=pval_col, effect_col=effect_col, stderr_col=stderr_col,
LD_reference=LD_reference, superpopulation=superpopulation)
before_after_plots(gene, susieDF, topVariants=topVariants)
before_after_consensus(gene, SumStats_sig, susieDF, max_SNPs=10)
# Create summary table for all genes
newEntry <- cbind(data.table(Gene=gene), subset(susieDF, PIP==max(PIP)) %>% as.data.table())
fineMapped_topSNPs <- rbind(fineMapped_topSNPs, newEntry)
cat("\n")
}
createDT(fineMapped_topSNPs, "Potential Causal SNPs Identified by susieR", scrollY = 200)
return(fineMapped_topSNPs)
}list[top_SNPs, SumStats_sig] <- import_sig_GWAS(
file_path = "Data/Parkinsons/Nalls2018_S2_SummaryStats.xlsx",
sheet="Data",
chrom_col = "CHR", position_col = "BP", snp_col="SNP",
pval_col="P, all studies", effect_col="Beta, all studies", gene_col="Nearest Gene",
caption= "Nalls et al. (2018) PD GWAS Summary Stats")
finemapped_PD <- finemap_geneList(top_SNPs, geneList=c("LRRK2","GBAP1","SNCA","VPS13C","GCH1"), #unique(top_SNPs$Gene)
filePath="Data/Parkinsons/META.PD.NALLS2014.PRS.tsv",
snp_col = "MarkerName", pval_col = "P.value")Extracting SNPs flanking lead SNP…
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-4890.90196469039” [1] “objective:-4890.86512255644” [1] “objective:-4890.86493793192” [1] “objective:-4890.86493700424”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
## Warning: Removed 1 rows containing missing values (geom_hline).
1 / 10 (10%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
Extracting SNPs flanking lead SNP…
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-1768.43682307539” [1] “objective:-1767.99485312494” [1] “objective:-1767.99377794051” [1] “objective:-1767.99377530418” [1] “objective:-1767.99377529783” [1] “objective:-1767.99377529781”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
## Warning: Removed 1 rows containing missing values (geom_hline).
0 / 10 (0%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
Extracting SNPs flanking lead SNP…
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-3587.57986639256” [1] “objective:-3584.84439532171” [1] “objective:-3584.84224690147” [1] “objective:-3584.84224525982” [1] “objective:-3584.84224525865” [1] “objective:-3584.84224525865”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
2 / 10 (20%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
Extracting SNPs flanking lead SNP…
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-3928.09335227999” [1] “objective:-3928.05140701022” [1] “objective:-3928.05134302287” [1] “objective:-3928.05134292561”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
## Warning: Removed 1 rows containing missing values (geom_hline).
1 / 10 (10%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
Extracting SNPs flanking lead SNP…
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-2952.66209303207” [1] “objective:-2952.62332682449” [1] “objective:-2952.62312855485” [1] “objective:-2952.62312754298”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
## Warning: Removed 1 rows containing missing values (geom_hline).
0 / 10 (0%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
list[top_SNPs, SumStats_sig] <- import_sig_GWAS(
file_path = "Data/Alzheimers/Posthuma/AD_target_SNP.xlsx",
sheet = 3,
chrom_col = "Chr", position_col = "bp", snp_col="SNP",
pval_col="P", effect_col="Z", gene_col="Gene",
caption= "Posthuma et al. (2018) AD GWAS Summary Stats")
finemapped_AD <- finemap_geneList(top_SNPs, geneList=c("CLU/PTK2B","APOE"), #unique(top_SNPs$Gene)
filePath="Data/Alzheimers/Posthuma/phase3.beta.se.hrc.txt",
effect_col = "BETA", stderr_col = "SE", position_col = "BP")## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
Fine mapping with SusieR… [1] “objective:-6159.30949844661” [1] “objective:-6158.9541676718” [1] “objective:-6158.95404180037” [1] “objective:-6158.95404175575”
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
1 / 10 (10%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
+ Creating LD matrix… LD Reference Panel = 1KG_Phase1 ped stats and snps stats have been set. ‘p’ has been set. ‘mu’ and ‘sigma’ have been set.
## Warning in log(susieDF$Effect): NaNs produced
## Warning in log(susieDF$Effect): NaNs produced
0 / 10 (0%) of SNPs of the SNPs in the summary stats were confirmed after fine-mapping.